Developer CD Series 2000 November: Tool Chest/Dev.CD Nov 00 TC Disk 1.toast/Sample Code/Devices and Hardware/Velocity Engine/VelEng Multiprecision/vfactor source/vgiants.c

C/C++ Source or Header | 2000-09-28 | 45.9 KB | 2,803 lines | [TEXT/CWIE]

/**************************************************************
 *
 *    vgiants.c
 *
 *    AltiVec-based giant integer package
 *  Derived from original giants.c project at
 *  Next Software, Inc.
 *
 *  This package is part of ongoing research in the
 *    Advanced Computation Group, Apple Computer.
 *
 *    c. 1999 Apple Computer, Inc.
 *    All Rights Reserved.
 *
 *
 *************************************************************/

/*
    Disclaimer:    IMPORTANT:  This Apple software is supplied to you by Apple Computer, Inc.
                ("Apple") in consideration of your agreement to the following terms, and your
                use, installation, modification or redistribution of this Apple software
                constitutes acceptance of these terms.  If you do not agree with these terms,
                please do not use, install, modify or redistribute this Apple software.

                In consideration of your agreement to abide by the following terms, and subject
                to these terms, Apple grants you a personal, non-exclusive license, under Apple’s
                copyrights in this original Apple software (the "Apple Software"), to use,
                reproduce, modify and redistribute the Apple Software, with or without
                modifications, in source and/or binary forms; provided that if you redistribute
                the Apple Software in its entirety and without modifications, you must retain
                this notice and the following text and disclaimers in all such redistributions of
                the Apple Software.  Neither the name, trademarks, service marks or logos of
                Apple Computer, Inc. may be used to endorse or promote products derived from the
                Apple Software without specific prior written permission from Apple.  Except as
                expressly stated in this notice, no other rights or licenses, express or implied,
                are granted by Apple herein, including but not limited to any patent rights that
                may be infringed by your derivative works or by other works in which the Apple
                Software may be incorporated.

                The Apple Software is provided by Apple on an "AS IS" basis.  APPLE MAKES NO
                WARRANTIES, EXPRESS OR IMPLIED, INCLUDING WITHOUT LIMITATION THE IMPLIED
                WARRANTIES OF NON-INFRINGEMENT, MERCHANTABILITY AND FITNESS FOR A PARTICULAR
                PURPOSE, REGARDING THE APPLE SOFTWARE OR ITS USE AND OPERATION ALONE OR IN
                COMBINATION WITH YOUR PRODUCTS.

                IN NO EVENT SHALL APPLE BE LIABLE FOR ANY SPECIAL, INDIRECT, INCIDENTAL OR
                CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
                GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
                ARISING IN ANY WAY OUT OF THE USE, REPRODUCTION, MODIFICATION AND/OR DISTRIBUTION
                OF THE APPLE SOFTWARE, HOWEVER CAUSED AND WHETHER UNDER THEORY OF CONTRACT, TORT
                (INCLUDING NEGLIGENCE), STRICT LIABILITY OR OTHERWISE, EVEN IF APPLE HAS BEEN
                ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/

#include "giantsdebug.h"
#if __MWERKS__
#include <altivec.h>
#endif
#include "vgiants.h"
#include "vecarith.h"
#include <Memory.h>
#include <Quickdraw.h>        // for random
#include <stddef.h>
#include <stdlib.h>
#include <math.h>
#include "giantsstack.h"
#include <segload.h>        // for ExitToShell

//////////////////////////////////////////////
// globals
//////////////////////////////////////////////

static giant    cur_recip = NULL;
static giant    cur_den = NULL;



//////////////////////////////////////////////

#if 0
#define BORROW_GIANT()    newgiantbits(DEFAULT_GIANT_BITS)
#define RETURN_GIANT(a)    disposegiant(a)
#else
#define BORROW_GIANT()    popg()
#define RETURN_GIANT(a)    pushg(1)
#endif

#define ADDVECS    addvecs
#define SUBVECS    subvecs

#define max(a, b) (a > b ? a : b)

#if GDEBUG

void ASSERT_GVALID(giant g)
{
    GAssert((g->giantSign == -1) || (g->giantSign == 1));

    if (((long)g->vectors & 0x0f) != 0) {
        DebugStr("\pvectors not 16-byte aligned.");
    }

    if (g->giantDigits > g->giantCapacity) {
        DebugStr("\pgiant digits overflow!");
    }

    if (g->giantDigits != 0) {
        vector unsigned long    veczero = (vector unsigned long)(0);

        if (vec_all_eq(g->vectors[g->giantDigits-1], veczero)) {
            DebugStr("\pMSV is zero");
        }
    }
}
#endif

static void FATAL_ERROR(long err)
{
#pragma unused(err)
    DebugStr("\pfatal error in vgiants");
    //ExitToShell();
}





gmatrix
newgmatrix(
    void
)
/* Allocates space for a gmatrix struct, but does not
 * allocate space for the giants. */
{
    return((gmatrix) malloc (4*sizeof(giant)));
}


GiantStructPtr    newgiantbits(unsigned long capacityBits)
{
    unsigned long    giantCapacity;
    unsigned long    allocSize;
    GiantStructPtr    pg;

    if (!capacityBits) capacityBits = 1;

    giantCapacity = (capacityBits + BITS_PER_VECTOR-1) / BITS_PER_VECTOR;
    allocSize = giantCapacity * BYTES_PER_VECTOR;
    pg = (GiantStructPtr)NewPtr(sizeof(GiantStruct));

    GAssert(pg != nil);

    pg->giantSign = 1;
    pg->giantDigits = 0;
    pg->giantCapacity = giantCapacity;

    pg->vectors = (vector unsigned long *)NewPtr(allocSize);

    GAssert(pg->vectors != nil);

    return pg;
}

GiantStructPtr    newgiant(unsigned long capacityShorts)
{
    return newgiantbits(capacityShorts << 4);
}


void disposegiant(GiantStructPtr g)
{
    ASSERT_GVALID(g);

    DisposePtr((Ptr)g->vectors);
    DisposePtr((Ptr)g);
}

static void reallocgdigits(giant g, unsigned long vectorcount, Boolean save)
{
    long                i;

    vector unsigned long *pNewVectors = (vector unsigned long *)NewPtr( vectorcount * BYTES_PER_VECTOR );

    GAssert(vectorcount >= g->giantDigits);

    GAssert(pNewVectors != nil);

    if (pNewVectors != nil) {
        if (save) {
            for (i=0; i<g->giantDigits; i++) {
                pNewVectors[i] = g->vectors[i];
            }
        }

        DisposePtr((Ptr)g->vectors);
        g->vectors = pNewVectors;
        g->giantCapacity = vectorcount;
    }
}



static void reallocg(giant g, unsigned long newbits, Boolean save)
{
    unsigned long    vectorcount =  (newbits + 127) / 128;

    reallocgdigits(g, vectorcount, save);
}

static void EnsureDigits(giant g, unsigned long newdigits, Boolean save)
{
    if (g->giantCapacity < newdigits) {
        reallocgdigits(g, newdigits, save);
    }
}

static void EnsureBits(giant g, unsigned long newbits, Boolean save)
{
    if (g->giantCapacity * BITS_PER_VECTOR < newbits) {
        reallocg(g, newbits, save);
    }
}

static void EnsureBytes(giant g, unsigned long newbytes, Boolean save)
{
    EnsureBits(g, newbytes << 3, save);
}

static void AddMSDLong(giant g, unsigned long digit)
{
    unsigned long *topvector;

    if (g->giantCapacity == g->giantDigits) {
        reallocgdigits(g, g->giantDigits+1, true);
    }

    topvector = (unsigned long*)&g->vectors[g->giantDigits];

    *topvector++ = 0;
    *topvector++ = 0;
    *topvector++ = 0;
    *topvector = digit;

    g->giantDigits++;
}

#if 1
static void justg(giant g)
{
    if (g->giantDigits) {
        unsigned long        counter = g->giantDigits;
        unsigned long        *pLongs = (unsigned long*)&g->vectors[g->giantDigits];    // point to one past last

        while (counter--) {
            if (*--pLongs) break;
            if (*--pLongs) break;
            if (*--pLongs) break;
            if (*--pLongs) break;

            g->giantDigits--;
        }

    }

}
#else
static void justg(giant g)
{
    if (g->giantDigits) {
        vector unsigned long vecZero = (vector unsigned long)(0);

        while (vec_all_eq(vecZero, g->vectors[g->giantDigits-1]) && g->giantDigits) {
            // do nothing
            g->giantDigits--;
        }
    }
}
#endif

static unsigned long gtoul(giant g)
{
    GAssert(bitlen(g) <= 32);

    if (!g->giantSign) return 0;

    return ((unsigned long *)g->vectors)[LONGS_PER_VECTOR-1];
}





int                cur_run = 0;
double *        sinCos=NULL;

void
init_sinCos(
    int         n
)
{
    int         j;
    double     e = TWOPI/n;

    if (n<=cur_run)
        return;
    cur_run = n;
    if (sinCos)
        free(sinCos);
    sinCos = (double *)malloc(sizeof(double)*(1+(n>>2)));
    for (j=0;j<=(n>>2);j++)
    {
        sinCos[j] = sin(e*j);
    }
}


double
s_sin(
    int     n
)
{
    int     seg = n/(cur_run>>2);

    switch (seg)
    {
        case 0: return(sinCos[n]);
        case 1: return(sinCos[(cur_run>>1)-n]);
        case 2: return(-sinCos[n-(cur_run>>1)]);
        case 3: return(-sinCos[cur_run-n]);
    }
    return 0;
}


double
s_cos(
    int     n
)
{
    int     quart = (cur_run>>2);

    if (n < quart)
        return(s_sin(n+quart));
    return(-s_sin(n-quart));
}




unsigned char GetNthByte(giant g, unsigned long pos)
{
    unsigned long    byteindex;

    //GAssert(pos < g->giantDigits * BYTES_PER_VECTOR);

    if (pos >= g->giantDigits * BYTES_PER_VECTOR) {
        return 0;
    }

    byteindex = pos & ~0xf;

    byteindex += BYTES_PER_VECTOR - (pos & 0xf)-1;

    return ((unsigned char*)g->vectors)[byteindex];

}

static unsigned short * GetNthShortPtr(giant g, unsigned long pos)
{
    unsigned long    shortindex;

    if (pos >= g->giantDigits * SHORTS_PER_VECTOR) {
        return 0;
    }

    shortindex = pos & ~0x7;

    shortindex += SHORTS_PER_VECTOR - (pos & 0x7)-1;

    return &((unsigned short*)g->vectors)[shortindex];

}


static unsigned long * GetNthLongPtr(giant g, unsigned long pos)
{
    unsigned long    longindex;

    if (pos >= g->giantDigits * LONGS_PER_VECTOR) {
        return 0;
    }

    longindex = pos & ~0x3;

    longindex += LONGS_PER_VECTOR - (pos & 0x3)-1;

    return &((unsigned long*)g->vectors)[longindex];

}

static unsigned long NumLongs(giant g)
{
    unsigned long result = 0;
    unsigned long *longPtr;
    unsigned long i;

    if (g->giantDigits) {
        result = g->giantDigits * LONGS_PER_VECTOR;

        longPtr = (unsigned long*)&g->vectors[g->giantDigits-1];

        for (i=0; i<LONGS_PER_VECTOR; i++) {
            if (*longPtr++ != 0)  break;
            result--;
        }
    }

    return result;
}


static unsigned long NumBytes(giant g)
{
    unsigned long result = 0;
    unsigned char *bytePtr;
    unsigned long i;

    if (g->giantDigits) {
        result = g->giantDigits * BYTES_PER_VECTOR;

        bytePtr = (unsigned char*)&g->vectors[g->giantDigits-1];

        for (i=0; i<BYTES_PER_VECTOR; i++) {
            if (*bytePtr++ != 0)  break;
            result--;
        }
    }

    return result;
}

static unsigned long BytesToDigits(unsigned long bytes)
{
    return (bytes+(BYTES_PER_VECTOR-1))/BYTES_PER_VECTOR;
}




static Boolean gshiftleftsubvector(unsigned long bits, giant g)
{
    vector unsigned char    temp;
    vector unsigned char    shiftFactor, negShiftFactor;
    vector unsigned char    vecZero;
    vector unsigned long    t1, t2;
    long                        i;
    Boolean                    shiftedOut = false;

    ASSERT_GVALID(g);

    if (bits) {
        vecZero = (vector unsigned char)(0);

        *((unsigned char*)&temp) = bits;

        shiftFactor = vec_splat(temp, 0);                // fill vector with shift amount
        negShiftFactor = vec_sub(vecZero, shiftFactor);

        // figure out if we need to expand
        t2 = vec_sro(g->vectors[g->giantDigits-1], negShiftFactor);
        t2 = vec_srl(t2, negShiftFactor);
        if (!vec_all_eq((vector unsigned long)vecZero, t2)) {
            shiftedOut = true;
            g->vectors[g->giantDigits] = t2;
        }

        for (i=g->giantDigits-1; i>0; i--) {
            t1 = vec_slo(g->vectors[i], shiftFactor);
            t1 = vec_sll(t1, shiftFactor);
            t2 = vec_sro(g->vectors[i-1], negShiftFactor);
            t2 = vec_srl(t2, negShiftFactor);
            g->vectors[i] = vec_or(t1, t2);
        }

        t1 = vec_slo(g->vectors[i], shiftFactor);
        g->vectors[i] = vec_sll(t1, shiftFactor);
    }

    return shiftedOut;
}

void         gshiftleft(unsigned long bits, giant g)
{
    unsigned long    bitlength;
    long            multiples;

    if (!bits) return;
    if (!g->giantDigits) return;

    bitlength = bitlen(g);

    if (bitlength + bits > (g->giantCapacity*BITS_PER_VECTOR)) {
        reallocg(g, bitlength + bits, true);
    }

    multiples = bits >> 7;

    if (multiples) {
        vector unsigned long    *pSrc;
        vector unsigned long    *pDst;
        unsigned long            digitCount = g->giantDigits;

        pSrc = &g->vectors[digitCount-1];
        pDst = &g->vectors[digitCount+multiples-1];

        g->giantDigits += multiples;

        // move up digits
        while (digitCount--) {
            *pDst-- = *pSrc--;
        }

        // fill remaining shifting with zeros
        while (multiples--) {
            *pDst-- = (vector unsigned long)(0);
        }
    }

    if (gshiftleftsubvector(bits & 0x7f, g)) {
        g->giantDigits++;
    }

    ASSERT_GVALID(g);
}


static void gshiftrightsubvector(unsigned long bits, giant g)
{
    vector unsigned char    temp;
    vector unsigned char    shiftFactor, negShiftFactor;
    vector unsigned char    vecZero;
    vector unsigned long    t1, t2;
    long                        i;

    if (bits) {
        vecZero = (vector unsigned char)(0);

        // set 0th element of shiftFactor to be bits
        *((unsigned char*)&temp) = bits;

        shiftFactor = vec_splat(temp, 0);    // fill vector with shift amount
        negShiftFactor = vec_sub(vecZero, shiftFactor);

        for (i=0; i<g->giantDigits-1; i++) {
            t1 = vec_sro(g->vectors[i], shiftFactor);
            t1 = vec_srl(t1, shiftFactor);
            t2 = vec_slo(g->vectors[i+1], negShiftFactor);
            t2 = vec_sll(t2, negShiftFactor);
            g->vectors[i] = vec_or(t1, t2);
        }

        t1 = vec_sro(g->vectors[i], shiftFactor);
        g->vectors[i] = vec_srl(t1, shiftFactor);
    }
}

void         gshiftright(unsigned long bits, giant g)
{
    long        multiples;

    if (!bits) return;
    if (!g->giantDigits) return;

    multiples = bits >> 7;

    if (multiples) {
        if (multiples >= g->giantDigits) {
            g->giantDigits = 0;
        } else {
            vector unsigned long    *pSrc;
            vector unsigned long    *pDst;
            unsigned long            newdigits = g->giantDigits-multiples;

            g->giantDigits = newdigits;

            pSrc = &g->vectors[multiples];
            pDst = g->vectors;

            while (newdigits--) {
                *pDst++ = *pSrc++;
            }
        }
    }


    if (g->giantDigits) {
        gshiftrightsubvector(bits& 0x7f, g);
    }

    justg(g);

    ASSERT_GVALID(g);
}




#if 1
unsigned long            bitlen(giant g)
{
    unsigned long            length;
    unsigned long             *longPtr;
    unsigned long            i;
    register unsigned long    longdigit;
    ASSERT_GVALID(g);

    if (!g->giantDigits) return 0;

    length = length = BITS_PER_VECTOR * g->giantDigits;

    longPtr = (unsigned long*)&g->vectors[g->giantDigits-1];

    for (i=0; i<4; i++) {
        longdigit = longPtr[i];

        if (!longdigit) {
            length -= 32;    // bits per long
        } else {
#if __MWERKS__
            register unsigned long leadingzeros=0;     // initialized to eliminate compiler warning
            asm {
                cntlzw    leadingzeros, longdigit
            }
            length -= leadingzeros;
#else
            while (!(longdigit & 0x80000000)) {
                longdigit = longdigit << 1;
                length -= 1;
            }

#endif
            break;
        }
    }

    return length;

}
#else
unsigned long            bitlen(giant g)
{
    unsigned long            length;
    unsigned long            count;
#if __MWERKS__
     vector unsigned long    vecones = (vector signed long)(-1);
#else
    vector unsigned long    vecones = (vector unsigned long)(-1);
#endif
    vector unsigned long    vecmask = (vector unsigned long)(0);
    vector unsigned long    veczero = (vector unsigned long)(0);
    vector unsigned long    digit;
    unsigned long            longindex;
    unsigned long            thelong;

    if (!g->giantDigits) return 0;

    count = 3;

    digit = g->vectors[g->giantDigits-1];

    longindex = 0;

    length = BITS_PER_VECTOR * g->giantDigits;

    vecmask = vec_sld(vecones, vecmask, 12);

    while (count-- && vec_all_eq(veczero, vec_and(vecmask, digit))) {
        vecmask = vec_sld(vecones, vecmask, 12);
        length -= 32;
        longindex++;
    }

    thelong = ((unsigned long*)(&g->vectors[g->giantDigits-1]))[longindex];

    while (!(thelong & 0x80000000)) {
        thelong = thelong << 1;
        length -= 1;
    }

    return length;
}
#endif

long    bitval(giant g, long pos)
{
#if 1
    return (GetNthByte(g, pos >> 3) & (1 << (pos & 0x07)));
#else
    unsigned char    mask;
    unsigned long    byteindex;

    byteindex = (pos >> 7) << 4;

    byteindex += BYTES_PER_VECTOR - ((pos & 0x7f) >> 3)-1;

    mask = 1 << (pos & 0x07);

    return ((unsigned char*)g->vectors)[byteindex] & mask;
#endif
}



void gtog(giant a, giant b)
{
    long i;

    ASSERT_GVALID(a);
    ASSERT_GVALID(b);

    if (b->giantCapacity < a->giantDigits) {
        reallocgdigits(b, a->giantDigits, false);
    }

    GAssert(b->giantCapacity >= a->giantDigits)


    for (i=0; i<a->giantDigits; i++) {
        b->vectors[i] = a->vectors[i];
    }

    b->giantDigits = a->giantDigits;
    b->giantSign = a->giantSign;
}


long CompareMagnitude(giant a, giant b)
{
    long                vectorindex, longindex;
    unsigned long         *aPtr, *bPtr;

    ASSERT_GVALID(a);
    ASSERT_GVALID(b);

    if (a->giantDigits > b->giantDigits) return 1;
    if (a->giantDigits < b->giantDigits) return -1;


    for (vectorindex=a->giantDigits-1; vectorindex >= 0; vectorindex--) {
        aPtr = (unsigned long*)&a->vectors[vectorindex];
        bPtr = (unsigned long*)&b->vectors[vectorindex];

        for (longindex = 0; longindex< LONGS_PER_VECTOR; longindex++) {
            if (*aPtr > *bPtr) return 1;
            if (*aPtr < *bPtr) return -1;
            aPtr++;
            bPtr++;
        }
    }

    return 0;
}



static void subg_normal(GiantStructPtr a, GiantStructPtr b, GiantStructPtr result)
{
    unsigned long endDigits = b->giantDigits;

    EnsureDigits(result, b->giantDigits, true);

    SUBVECS(a->vectors, a->giantDigits, b->vectors, b->giantDigits, result->vectors);

    result->giantDigits = endDigits;

    justg(result);
}


void addg(giant a, giant b)
{
    ASSERT_GVALID(a);
    ASSERT_GVALID(b);

    if (isZero(a)) {
        return;
    }

    if (isZero(b)) {
        gtog(a, b);
        return;
    }

    if (a->giantSign == b->giantSign) {
        unsigned long    newdigits;

        newdigits = max(a->giantDigits, b->giantDigits)+1;
        EnsureDigits(b, newdigits, true);
        if (!ADDVECS(a->vectors, a->giantDigits, b->vectors, b->giantDigits, b->vectors)) {
            newdigits--;
        }
        b->giantDigits = newdigits;

    } else {
        if (b->giantSign > 0) {
            // a is negative, b is positive
            if (CompareMagnitude(b, a) >= 0) {
                subg_normal(a, b, b);
            } else {
                // a neg, b pos, b<a
                subg_normal(b, a, b);
                negg(b);
            }
        } else {
            // a is positive, b is negative
            if (CompareMagnitude(b, a) < 0) {
                subg_normal(b, a, b);
                negg(b);
            } else {
                // a neg, b pos, b>a
                subg_normal(a, b, b);
            }
        }
    }
    ASSERT_GVALID(a);
    ASSERT_GVALID(b);
}



// result = b-a
void subg(giant a, giant b)
{
    ASSERT_GVALID(a);
    ASSERT_GVALID(b);

    if (isZero(a)) {
        return;
    }

    if (isZero(b)) {
        gtog(a, b);
        negg(b);
        return;
    }

    if (a->giantSign != b->giantSign) {
        unsigned long    newdigits;

        newdigits = max(a->giantDigits, b->giantDigits)+1;
        EnsureDigits(b, newdigits, true);
        if (!ADDVECS(a->vectors, a->giantDigits, b->vectors, b->giantDigits, b->vectors)) {
            newdigits--;
        }
        b->giantDigits = newdigits;
    } else {
        if (b->giantSign > 0) {
            // a is negative, b is positive
            if (CompareMagnitude(b, a) >= 0) {
                subg_normal(a, b, b);
            } else {
                // a neg, b pos, b<a
                subg_normal(b, a, b);
                negg(b);
            }
        } else {
            // a is positive, b is negative
            if (CompareMagnitude(b, a) < 0) {
                subg_normal(b, a, b);
                negg(b);
            } else {
                // a neg, b pos, b>a
                subg_normal(a, b, b);
            }
        }
    }

    ASSERT_GVALID(a);
    ASSERT_GVALID(b);
}





void grammarmulg(giant a, giant b)
{
    unsigned long    resultdigits = a->giantDigits + b->giantDigits;

    if (isZero(b)) return;

    if (isZero(a)) {
        gtog(a, b);
        return;
    }

    if (a->giantDigits == 1) {
        if (NumLongs(a) == 1) {

            smulg(gtoul(a), b);
            b->giantSign *= a->giantSign;

        } else {
            EnsureDigits(b, b->giantDigits+1, true);

            VecMult1(    a->vectors,
                        b->giantDigits,
                        b->vectors,
                        b->vectors);

            b->giantDigits++;

            b->giantSign *= a->giantSign;

            justg(b);
        }
    } else  if (b->giantDigits == 1) {
        if (NumLongs(b) == 1) {
            unsigned long     blong = gtoul(b);
            long             bsign = b->giantSign;

            gtog(a, b);
            smulg(blong, b);
            b->giantSign *= bsign;

        } else {
            EnsureDigits(b, a->giantDigits+1, true);

            VecMult1(    b->vectors,
                        a->giantDigits,
                        a->vectors,
                        b->vectors);

            b->giantDigits = a->giantDigits+1;

            b->giantSign *= a->giantSign;

            justg(b);
        }
    } else {
        GiantStructPtr     temp = BORROW_GIANT();

        if (a->giantDigits == 2 && b->giantDigits == 2) {
            EnsureDigits(temp, 4, false);

            mult2by2(    &a->vectors[1],
                        &a->vectors[0],
                        &b->vectors[1],
                        &b->vectors[0],
                        temp->vectors);

            temp->giantDigits = 4;

            temp->giantSign = a->giantSign*b->giantSign;

            justg(temp);

            gtog(temp, b);

            RETURN_GIANT(temp);

        } else {
            EnsureDigits(temp, resultdigits, false);

            VecMult(a->vectors,
                    a->giantDigits,
                    b->vectors,
                    b->giantDigits,
                    temp->vectors);

            temp->giantDigits = resultdigits;

            temp->giantSign = a->giantSign*b->giantSign;
            justg(temp);

            gtog(temp, b);

            RETURN_GIANT(temp);
        }
    }

    ASSERT_GVALID(b);

}

void
karatmulg(
        giant   x,
        giant   y)
/* y becomes x*y. */
{
        int     s = x->giantDigits, t = y->giantDigits, w, bits,
                sg = gsign(x)*gsign(y);
        giant   a, b, c, e, f;

        if((s <= KARAT_BREAK) || (t <= KARAT_BREAK)) {
                grammarmulg(x,y);
                return;
        }
        w = (s + t + 2)/4; bits = 16*w;

        a = BORROW_GIANT(); b = BORROW_GIANT(); c = BORROW_GIANT(); e = BORROW_GIANT(); f = BORROW_GIANT();
        gtog(x,a); if (w <= s) {a->giantDigits = w; justg(a);}
        absg(a);
        gtog(x,b); absg(b);
        gshiftright(bits, b);

        gtog(y,c); if (w <= t) {c->giantDigits = w; justg(c);}
        absg(c);
        absg(y);
        gshiftright(bits,y);
        gtog(a,e); addg(b,e);
        gtog(c,f); addg(y,f);
        karatmulg(e,f);  /* f := (a + b)(c + d). */
        karatmulg(c,a); /* a := a c. */
        karatmulg(b,y); /* d := b d. */
        subg(a,f); subg(y,f);
        gshiftleft(bits, y);
        addg(f, y);
        gshiftleft(bits, y);

        addg(a, y);
        setsign(y, sg);

        RETURN_GIANT(f);
        RETURN_GIANT(e);
        RETURN_GIANT(c);
        RETURN_GIANT(b);
        RETURN_GIANT(a);

        return;
}

void mulg(giant a, giant b)
{
    if((a->giantDigits <= KARAT_BREAK) || (b->giantDigits <= KARAT_BREAK)) {
        grammarmulg(a,b);
    } else {
        karatmulg(a, b);
    }
}


#if 1
/* g += i, with i non-negative. */
void         iaddg(unsigned long i,giant g)
{
    if (isZero(g)) {
        itog(i, g);
    } else {
        ASSERT_GVALID(g);

        if (AddULongToVecs(g->vectors, g->giantDigits, i)) {
            AddMSDLong(g, 1);
        }

        justg(g);

        ASSERT_GVALID(g);
    }
}
#else
/* g += i, with i non-negative. */
void         iaddg(unsigned long i,giant g)
{
    giant        temp=BORROW_GIANT();

    ASSERT_GVALID(g);

    itog(i, temp);
    addg(temp, g);

    ASSERT_GVALID(g);

    RETURN_GIANT(temp);
}
#endif



/* Integer <-> giant. */
void itog(long n, giant g)
{
    unsigned long    absn = abs(n);

    if (absn) {
        unsigned long    *pLong = (unsigned long *)g->vectors;

        *pLong++ = 0;
        *pLong++ = 0;
        *pLong++ = 0;
        *pLong = absn;

        g->giantSign = (n < 0 ? -1 : 1);
        g->giantDigits = 1;

    } else {
        g->giantDigits = 0;
    }

    ASSERT_GVALID(g);
}


#if 1
void         squareg(giant g)
{
    GiantStructPtr    temp;
    unsigned long     digits;

    if (isZero(g)) return;

    temp = BORROW_GIANT();

    digits = g->giantDigits * 2;

    EnsureDigits(temp, digits, false);

    VecSquare(g->vectors,
            g->giantDigits,
            temp->vectors);

    temp->giantDigits = digits;

    temp->giantSign = 1;

    justg(temp);

    gtog(temp, g);

    RETURN_GIANT(temp);

    ASSERT_GVALID(g);
}

#else
/* g *= g. */
void         squareg(giant g)
{
    ASSERT_GVALID(g);
    mulg(g, g);

    ASSERT_GVALID(g);

}
#endif

void         randomg(giant g, unsigned long numvecs)
{
    long i;

    if (numvecs > g->giantCapacity) numvecs = g->giantCapacity;

    for (i=0; i<numvecs*SHORTS_PER_VECTOR; i++) {
        ((unsigned short*)g->vectors)[i] = Random();
    }

    g->giantSign = (Random() & 1) ? -1 : 1;
    g->giantDigits = numvecs;

    justg(g);

    ASSERT_GVALID(g);
}




void
modg_via_recip(
    giant     d,
    giant     r,
    giant     n
)
/* This is the fastest mod of the present collection.
 * n := n % d, where r is the precalculated
 * steady-state reciprocal of d. */

{
    int        s = (bitlen(r)-1), sign = gsign(n);
    giant     tmp, tmp2;

    if (isZero(d) || (gsign(d) < 0))
    {
        FATAL_ERROR(SIGN);
    }

    tmp = BORROW_GIANT();
    tmp2 = BORROW_GIANT();

    if (isone(d)) {
        itog(0, n);
        goto done;
    }

    absg(n);
    while (1)
    {
        gtog(n, tmp); gshiftright(s-1, tmp);
        mulg(r, tmp);
        gshiftright(s+1, tmp);
        mulg(d, tmp);
        subg(tmp, n);
        if (gcompg(n,d) >= 0)
            subg(d,n);
        if (gcompg(n,d) < 0)
            break;
    }
    if (sign >= 0)
        goto done;
    if (isZero(n))
        goto done;
    negg(n);
    addg(d,n);
done:
    RETURN_GIANT(tmp);
    RETURN_GIANT(tmp2);
    return;
}

void
make_recip(
    giant     d,
    giant     r
)
/* r becomes the steady-state reciprocal
 * 2^(2b)/d, where b = bit-length of d-1. */
{
    int        b;
    giant     tmp, tmp2;

    if (isZero(d) || (gsign(d) < 0))
    {
        FATAL_ERROR(SIGN);
    }
    tmp = BORROW_GIANT();
    tmp2 = BORROW_GIANT();
    itog(1, r);
    subg(r, d);
    b = bitlen(d);
    addg(r, d);
    gshiftleft(b, r);
    gtog(r, tmp2);
    while (1)
    {
        gtog(r, tmp);
        squareg(tmp);
        gshiftright(b, tmp);
        mulg(d, tmp);
        gshiftright(b, tmp);
        addg(r, r);
        subg(tmp, r);
        if (gcompg(r, tmp2) <= 0)
            break;
        gtog(r, tmp2);
    }
    itog(1, tmp);
    gshiftleft(2*b, tmp);
    gtog(r, tmp2);
    mulg(d, tmp2);
    subg(tmp2, tmp);
    itog(1, tmp2);
    while (gsign(tmp) < 0)
    {
        subg(tmp2, r);
        addg(d, tmp);
    }
    RETURN_GIANT(tmp);
    RETURN_GIANT(tmp2);
}


void
modg(
    giant     d,
    giant     n
)
/* n becomes n%d. n is arbitrary, but the denominator d must be positive! */
{
    if (cur_recip == NULL) {
        cur_recip = newgiantbits(1);
        cur_den = newgiantbits(1);
        gtog(d, cur_den);
        make_recip(d, cur_recip);
    } else if (gcompg(d, cur_den)) {
        gtog(d, cur_den);
        make_recip(d, cur_recip);
    }
    modg_via_recip(d, cur_recip, n);
}

void
divg_via_recip(
    giant     d,
    giant     r,
    giant     n
)
/* n := n/d, where r is the precalculated
 * steady-state reciprocal of d. */
{
    int     s = 2*(bitlen(r)-1), sign = gsign(n);
    giant     tmp, tmp2;

    if (isZero(d) || (gsign(d) < 0))
    {
        FATAL_ERROR(SIGN);
    }

    tmp = BORROW_GIANT();
    tmp2 = BORROW_GIANT();

    if (isone(d)) {
        // dividing by one!
        goto done;
    }

    absg(n);
    itog(0, tmp2);
    while (1)
    {
        gtog(n, tmp);
        mulg(r, tmp);
        gshiftright(s, tmp);
        addg(tmp, tmp2);
        mulg(d, tmp);
        subg(tmp, n);
        if (gcompg(n,d) >= 0)
        {
            subg(d,n);
            iaddg(1, tmp2);
        }
        if (gcompg(n,d) < 0)
            break;
    }
    gtog(tmp2, n);
    setsign(n, gsign(n) * sign);
done:
    RETURN_GIANT(tmp);
    RETURN_GIANT(tmp2);
}



void
divg(
    giant     d,
    giant     n
)
/* n becomes n/d. n is arbitrary, but the denominator d must be positive! */
{
    if (cur_recip == NULL) {
        cur_recip = newgiantbits(1);
        cur_den = newgiantbits(1);
        gtog(d, cur_den);
        make_recip(d, cur_recip);
    } else if (gcompg(d, cur_den)) {
        gtog(d, cur_den);
        make_recip(d, cur_recip);
    }
    divg_via_recip(d, cur_recip, n);
}

static long
radixdiv(
    long                base,
    long                divisor,
    giant            thegiant
)
/* Divides giant of arbitrary base by divisor.
 * Returns remainder. Used by idivg and gread. */
{
    long                j = (thegiant->giantDigits*SHORTS_PER_VECTOR)-1;
    unsigned int         num,rem=0;
    unsigned short        *digitpointer;

    if (divisor == 0)
    {
        FATAL_ERROR(DIVIDEBYZERO);
    }

    while (j>=0)
    {
        digitpointer = GetNthShortPtr(thegiant, j);
        num=rem*base + *digitpointer;
        *digitpointer = (unsigned short)(num/divisor);

        rem = num % divisor;
        --j;
    }

    if (divisor < 0) negg(thegiant);

    justg(thegiant);

    return(rem);
}

#if 1
long
idivg(
    long    divisor,
    giant     theg
)
{
    /* theg becomes theg/divisor. Returns remainder. */
    int     n;
    int     base = 1<<(8*sizeof(short));

    n = radixdiv(base,divisor,theg);
    return(n);
}
#else
long
idivg(
    long        divisor,
    giant         theg
)
{
    giant    originalg = BORROW_GIANT();
    giant     divisorg = BORROW_GIANT();
    long    remainder;

    gtog(theg, originalg);

    itog(divisor, divisorg);

    divg(divisorg, theg);

    mulg(theg, divisorg);

    subg(divisorg, originalg);

    remainder = gtoi(originalg);

    RETURN_GIANT(divisorg);
    RETURN_GIANT(originalg);

    return remainder;
}
#endif





static void
gswap(
    giant     *p,
    giant     *q
)
{
    giant     t;

    t = *p;
    *p = *q;
    *q = t;
}


static void
fix(
    giant     *p,
    giant     *q
)
/* Insure that x > y >= 0. */
{
    if( gsign(*p) < 0 )
        negg(*p);
    if( gsign(*q) < 0 )
        negg(*q);
    if( gcompg(*p,*q) < 0 )
        gswap(p,q);
}


unsigned long    numtrailzeros(giant g)
{
    unsigned long trailingbits = 0;

    ASSERT_GVALID(g);

    if (g->giantDigits != 0) {
        vector unsigned long    veczero = (vector unsigned long)(0);
        long                    vecindex;
        long                    byteindex;
        long                    bitindex;
        unsigned char            lastbyte;

        for (vecindex = 0; vecindex < g->giantDigits; vecindex++) {
            if (!vec_all_eq(veczero, g->vectors[vecindex])) {
                break;
            }

            trailingbits += BITS_PER_VECTOR;
        }

        // now go by bytes;
        for (byteindex = 0; byteindex < BYTES_PER_VECTOR; byteindex++) {
            lastbyte = GetNthByte(g, byteindex + (vecindex * BYTES_PER_VECTOR));
            if (lastbyte != 0) {
                break;
            }

            trailingbits += 8;
        }


        for (bitindex=0; bitindex < 8; bitindex++) {
            if (lastbyte & 1) {
                break;
            }

            lastbyte = lastbyte >> 1;
            trailingbits++;
        }
    }

    return trailingbits;
}

static void
dotproduct(
    giant     a,
    giant     b,
    giant     c,
    giant     d
)
/* Replace last argument with the dot product of two 2-vectors. */
{
    giant s4 = BORROW_GIANT();

    gtog(c,s4);
    mulg(a, s4);
    mulg(b,d);
    addg(s4,d);

    RETURN_GIANT(s4);
}


static   void
mulvM(
    gmatrix     A,
    giant         x,
    giant         y
)
/* Multiply vector by Matrix; changes x,y. */
{
    giant s0 = BORROW_GIANT();
    giant s1 = BORROW_GIANT();

    gtog(A->ur,s0);
    gtog( A->lr,s1);
    dotproduct(x,y,A->ul,s0);
    dotproduct(x,y,A->ll,s1);
    gtog(s0,x);
    gtog(s1,y);

    RETURN_GIANT(s0);
    RETURN_GIANT(s1);
}


static void
bgcdg(
    giant     a,
    giant     b
)
/* Binary form of the gcd. b becomes the gcd of a,b. */
{
    int        k = isZero(b), m = isZero(a);
    giant     u, v, t;

    if (k || m)
    {
        if (m)
        {
            if (k)
                itog(1,b);
            return;
        }
        if (k)
        {
            if (m)
                itog(1,b);
            else
                gtog(a,b);
            return;
        }
    }

    u = BORROW_GIANT();
    v = BORROW_GIANT();
    t = BORROW_GIANT();

    /* Now neither a nor b is zero. */
    gtog(a, u);
    absg(u);
    gtog(b, v);
    absg(v);
    k = numtrailzeros(u);
    m = numtrailzeros(v);
    if (k>m)
        k = m;
    gshiftright(k,u);
    gshiftright(k,v);
    if (isOdd(u)) {
        gtog(v, t);
        negg(t);
    }
    else
    {
        gtog(u,t);
    }

    while (!isZero(t))
    {
        m = numtrailzeros(t);
        gshiftright(m, t);
        if (gsign(t) > 0)
        {
            gtog(t, u);
            subg(v,t);
        }
        else
        {
            gtog(t, v);
            negg(v);
            addg(u,t);
        }
    }
    gtog(u,b);
    gshiftleft(k, b);
    RETURN_GIANT(v);
    RETURN_GIANT(u);
    RETURN_GIANT(t);
}

static void
shgcd(
    register int    x,
    register int     y,
    gmatrix         A
)
/*
 * Do a half gcd on the integers a and b, putting the result in A
 * It is fairly easy to use the 2 by 2 matrix description of the
 * extended Euclidean algorithm to prove that the quantity q*t
 * never overflows.
 */
{
    register int     q, t, start = x;
    int             Aul = 1, Aur = 0, All = 0, Alr = 1;

    while    (y != 0 && y > start/y)
    {
        q = x/y;
        t = y;
        y = x%y;
        x = t;
        t = All;
        All = Aul-q*t;
        Aul = t;
        t = Alr;
        Alr = Aur-q*t;
        Aur = t;
    }
    itog(Aul,A->ul);
    itog(Aur,A->ur);
    itog(All,A->ll);
    itog(Alr,A->lr);
}

static void
punch(
    giant         q,
    gmatrix     A
)
/* Multiply the matrix A on the left by [0,1,1,-q]. */
{
    giant s0 = BORROW_GIANT();

    gtog(A->ll,s0);
    mulg(q,A->ll);
    gswap(&A->ul,&A->ll);
    subg(A->ul,A->ll);
    gtog(s0,A->ul);
    gtog(A->lr,s0);
    mulg(q,A->lr);
    gswap(&A->ur,&A->lr);
    subg(A->ur,A->lr);
    gtog(s0,A->ur);

    RETURN_GIANT(s0);
}


static void
onestep(
    giant         x,
    giant         y,
    gmatrix     A
)
/* Do one step of the euclidean algorithm and modify
 * the matrix A accordingly. */
{
    giant s4 = BORROW_GIANT();

    gtog(x,s4);
    gtog(y,x);
    gtog(s4,y);
    divg(x,s4);
    punch(s4,A);
    mulg(x,s4);
    subg(s4,y);

    RETURN_GIANT(s4);
}

static void
mulmM(
    gmatrix     A,
    gmatrix     B
)
/* Multiply matrix by Matrix; changes second matrix. */
{
    giant s0 = BORROW_GIANT();
    giant s1 = BORROW_GIANT();
    giant s2 = BORROW_GIANT();
    giant s3 = BORROW_GIANT();

    gtog(B->ul,s0);
    gtog(B->ur,s1);
    gtog(B->ll,s2);
    gtog(B->lr,s3);

    dotproduct(A->ur,A->ul,B->ll,s0);
    dotproduct(A->ur,A->ul,B->lr,s1);
    dotproduct(A->ll,A->lr,B->ul,s2);
    dotproduct(A->ll,A->lr,B->ur,s3);

    gtog(s0,B->ul);
    gtog(s1,B->ur);
    gtog(s2,B->ll);
    gtog(s3,B->lr);

    RETURN_GIANT(s0);
    RETURN_GIANT(s1);
    RETURN_GIANT(s2);
    RETURN_GIANT(s3);
}



static void
hgcd(
    int         n,
    giant         xx,
    giant         yy,
    gmatrix     A
)
/* hgcd(n,x,y,A) chops n bits off x and y and computes th
 * 2 by 2 matrix A such that A[x y] is the pair of terms
 * in the remainder sequence starting with x,y that is
 * half the size of x. Note that the argument A is modified
 * but that the arguments xx and yy are left unchanged.
 */
{
    giant         x, y;

    if (isZero(yy))
        return;

    x = BORROW_GIANT();
    y = BORROW_GIANT();
    gtog(xx,x);
    gtog(yy,y);
    gshiftright(n,x);
    gshiftright(n,y);
    if (bitlen(x) <= INTLIMIT )
    {
        shgcd(gtoi(x),gtoi(y),A);
    }
    else
    {
        GMatrixStruct     B;
        int             m = bitlen(x)/2;

        hgcd(m,x,y,A);
        mulvM(A,x,y);
        if (gsign(x) < 0)
        {
            negg(x); negg(A->ul); negg(A->ur);
        }
        if (gsign(y) < 0)
        {
            negg(y); negg(A->ll); negg(A->lr);
        }
        if (gcompg(x,y) < 0)
        {
            gswap(&x,&y);
            gswap(&A->ul,&A->ll);
            gswap(&A->ur,&A->lr);
        }
        if (!isZero(y))
        {
            onestep(x,y,A);
            m /= 2;
            B.ul = BORROW_GIANT();
            B.ur = BORROW_GIANT();
            B.ll = BORROW_GIANT();
            B.lr = BORROW_GIANT();
            itog(1,B.ul);
            itog(0,B.ur);
            itog(0,B.ll);
            itog(1,B.lr);
            hgcd(m,x,y,&B);
            mulmM(&B,A);

            RETURN_GIANT(B.ul);
            RETURN_GIANT(B.ur);
            RETURN_GIANT(B.ll);
            RETURN_GIANT(B.lr);
        }
    }
    RETURN_GIANT(x);
    RETURN_GIANT(y);
}




static void
ggcd(
    giant     xx,
    giant     yy
)
/* A giant gcd.  Modifies its arguments. */
{
    giant             x = BORROW_GIANT();
    giant            y = BORROW_GIANT();
    GMatrixStruct     A;

    gtog(xx,x); gtog(yy,y);
    for(;;)
    {
        fix(&x,&y);
        if (bitlen(y) <= GCDLIMIT )
            break;
        A.ul = BORROW_GIANT();
        A.ur = BORROW_GIANT();
        A.ll = BORROW_GIANT();
        A.lr = BORROW_GIANT();
        itog(1,A.ul);
        itog(0,A.ur);
        itog(0,A.ll);
        itog(1,A.lr);
        hgcd(0,x,y,&A);
        mulvM(&A,x,y);
        RETURN_GIANT(A.ul);
        RETURN_GIANT(A.ur);
        RETURN_GIANT(A.ll);
        RETURN_GIANT(A.lr);
        fix(&x,&y);
        if (bitlen(y) <= GCDLIMIT )
            break;
        modg(y,x);
        gswap(&x,&y);
    }
    bgcdg(x,y);
    gtog(y,yy);
    RETURN_GIANT(x);
    RETURN_GIANT(y);
}

void
gcdg(
    giant        x,
    giant        y
)
{
    if (bitlen(y)<= GCDLIMIT)
        bgcdg(x,y);
    else
        ggcd(x,y);
}


void
cgcdg(
    giant     a,
    giant     v
)
/* Classical Euclid GCD. v becomes gcd(a, v). */
{
    giant     u, r;

    absg(v);
    if (isZero(a))
        return;

    u = BORROW_GIANT();
    r = BORROW_GIANT();

    gtog(a, u);
    absg(u);
    while (!isZero(v))
    {
        gtog(u, r);
        modg(v, r);
        gtog(v, u);
        gtog(r, v);
    }
    gtog(u,v);
    RETURN_GIANT(u);
    RETURN_GIANT(r);
}




static void
columnwrite(
    FILE     *filepointer,
    short     *column,
    char     *format,
    short     arg,
    int     newlines
)
/* Used by gwriteln. */
{
    char     outstring[10];
    short     i;

    sprintf(outstring,format,arg);
    for (i=0; outstring[i]!=0; ++i)
    {
        if (newlines)
        {
            if (*column >= COLUMNWIDTH)
            {
                fputs("\\\n",filepointer);
                *column = 0;
            }
        }
        fputc(outstring[i],filepointer);
        ++*column;
    }
}


void
gwrite(
    giant            thegiant,
    FILE            *filepointer,
    int                newlines
)
/* Outputs thegiant to filepointer. Output is terminated by a newline. */
{
    short            column;
    static giant    garbagegiant = NULL;
    unsigned long    *pLongStorage;
    unsigned long    numlongs = (thegiant->giantDigits * 40);
    unsigned long    *numpointer;
    Boolean            firstdigits=true;

    pLongStorage = (unsigned long *)NewPtr(numlongs*sizeof(long));

    if (garbagegiant == NULL)
        garbagegiant = newgiantbits(1);

    if (isZero(thegiant))
    {
        fputs("0",filepointer);
    }
    else
    {
        numpointer = pLongStorage;
        gtog(thegiant,garbagegiant);
        absg(garbagegiant);

        do
        {
            *numpointer = (unsigned short)idivg(10000,garbagegiant);
            ++numpointer;
        }     while (!isZero(garbagegiant));

        column = 0;
        numpointer--;

        if (gsign(thegiant)<0) {
            columnwrite(filepointer,&column,"-",0, newlines);
        }

        do {
            if (firstdigits) {
                columnwrite(filepointer,&column,"%d",*numpointer--,newlines);
                firstdigits = false;
            } else {
                columnwrite(filepointer,&column,"%04d",*numpointer--,newlines);
            }
        } while (numpointer >= pLongStorage);

    }

    DisposePtr((Ptr)pLongStorage);
}

void
gwriteln(
    giant        theg,
    FILE        *filepointer
)
{
    gwrite(theg, filepointer, 1);
    fputc('\n',filepointer);
}


void
gread(
    giant             theg,
    FILE             *filepointer
)
{
    char             currentchar;
    int             isneg,size,backslash=false,numdigits=0;
    static giant    basetenthousand = NULL;
    static char        *inbuf = NULL;
    static giant    currentmult;

    if (basetenthousand == NULL)
        basetenthousand = newgiantbits(1);
    if (currentmult == NULL)
        currentmult = newgiantbits(1);
    if (inbuf == NULL)
        inbuf = (char*)malloc(MAX_DIGITS);

    itog(1, currentmult);
    itog(0, theg);

    currentchar = (char)fgetc(filepointer);
    if (currentchar=='-')
    {
        isneg=true;
    }
    else
    {
        isneg=false;
        if (currentchar!='+')
            ungetc(currentchar,filepointer);
    }

    do
    {
        currentchar = (char)fgetc(filepointer);
        if ((currentchar>='0') && (currentchar<='9'))
        {
            inbuf[numdigits]=currentchar;
            if(++numdigits==MAX_DIGITS)
                break;
            backslash=false;
        }
        else
        {
            if (currentchar=='\\')
                backslash=true;
        }
    } while(((currentchar!=' ') && (currentchar!='\n') &&
                (currentchar!='\t')) || (backslash) );
    if (numdigits)
    {
        size = 0;
        do
        {
            inbuf[numdigits] = 0;
            numdigits-=4;
            if (numdigits<0)
                numdigits=0;

            itog((unsigned short)strtol(&inbuf[numdigits],NULL,10), basetenthousand);
            mulg(currentmult, basetenthousand);
            addg(basetenthousand, theg);
            itog(10000, basetenthousand);
            mulg(basetenthousand, currentmult);

        } while (numdigits>0);


        if (isneg) {
            setsign(theg, -1);
        } else {
            setsign(theg, 1);
        }
    }
}

void         gdumphex(giant g)
{
    long            i;
    unsigned char    b;

    if (g->giantDigits) {
        if (g->giantSign == -1) {
            printf("-");
        }

        for (i=BYTES_PER_VECTOR*(g->giantDigits)-1; i>=0; i--) {
            b = GetNthByte(g, i);
            if (b == 0) {
                printf("00");
            } else if (b<0x10) {
                printf("0%x", b);
            } else {
                printf("%x", b);
            }
        }
    } else {
        printf("0");
    }
}



void
bdivg(
    giant        v,
    giant        u
)
/* u becomes greatest power of two not exceeding u/v. */
{
    int         diff = bitlen(u) - bitlen(v);
    giant        scratch7;

    if (diff<0)
    {
        itog(0,u);
        return;
    }
    scratch7 = BORROW_GIANT();
    gtog(v, scratch7);
    gshiftleft(diff,scratch7);
    if (gcompg(u,scratch7) < 0)
        diff--;
    if (diff<0)
    {
        itog(0,u);
        RETURN_GIANT(scratch7);
        return;
    }
    itog(1,u);
    gshiftleft(diff,u);

    RETURN_GIANT(scratch7);
}



static giant    u0=NULL, u1=NULL, v0=NULL, v1=NULL;


static int
binvaux(
    giant     p,
    giant     x
)
/* Binary inverse method. Returns zero if no inverse exists,
 * in which case x becomes GCD(x,p). */
{

    giant scratch7;

    if (isone(x))
        return(1);
    if (!v0)
        v0 = newgiantbits(1);
    if (!v1)
        v1 = newgiantbits(1);
    if (!u1)
        u1 = newgiantbits(1);
    if (!u0)
        u0 = newgiantbits(1);
    itog(1, v0);
    gtog(x, v1);
    itog(0,x);
    gtog(p, u1);

    scratch7 = BORROW_GIANT();
    while(!isZero(v1))
    {
        gtog(u1, u0);
        bdivg(v1, u0);
        gtog(x, scratch7);
        gtog(v0, x);
        mulg(u0, v0);
        subg(v0,scratch7);
        gtog(scratch7, v0);

        gtog(u1, scratch7);
        gtog(v1, u1);
        mulg(u0, v1);
        subg(v1,scratch7);
        gtog(scratch7, v1);
    }

    RETURN_GIANT(scratch7);

    if (!isone(u1))
    {
        gtog(u1,x);
        if(gsign(x)<0) addg(p, x);
            return(0);
    }
    if(gsign(x)<0)
        addg(p, x);
    return(1);
}

int
binvg(
    giant     p,
    giant     x
)
{
    modg(p, x);
    return(binvaux(p,x));
}

static int
invaux(
    giant     p,
    giant     x
)
/* Returns zero if no inverse exists, in which case x becomes
 * GCD(x,p). */
{

    giant scratch7;

    if (isone(x))
        return(1);
    if (!v0)
        v0 = newgiantbits(1);
    if (!v1)
        v1 = newgiantbits(1);
    if (!u1)
        u1 = newgiantbits(1);
    if (!u0)
        u0 = newgiantbits(1);
    itog(1,u1);
    gtog(p, v0);
    gtog(x, v1);
    itog(0,x);

    scratch7 = BORROW_GIANT();
    while (!isZero(v1))
    {
        gtog(v0, u0);
        divg(v1, u0);
        gtog(u0, scratch7);
        mulg(v1, scratch7);
        subg(v0, scratch7);
        negg(scratch7);
        gtog(v1, v0);
        gtog(scratch7, v1);
        gtog(u1, scratch7);
        mulg(u0, scratch7);
        subg(x, scratch7);
        negg(scratch7);
        gtog(u1,x);
        gtog(scratch7, u1);
    }

    RETURN_GIANT(scratch7);

    if (!isone(v0))
    {
        gtog(v0,x);
        return(0);
    }
    if(gsign(x)<0)
        addg(p, x);
    return(1);
}


int
invg(
    giant     p,
    giant     x
)
{
    modg(p, x);
    return(invaux(p,x));
}

void
mersennemod (
    int n,
    giant g
)
/* g := g (mod 2^n - 1) */
{
    int the_sign, s;
    giant scratch3 = BORROW_GIANT();
    giant scratch4 = BORROW_GIANT();

    if ((the_sign = gsign(g)) < 0) absg(g);
    while (bitlen(g) > n) {
        gtog(g,scratch3);
        gshiftright(n,scratch3);
        addg(scratch3,g);
        gshiftleft(n,scratch3);
        subg(scratch3,g);
    }
    if(!isZero(g)) {
        if ((s = gsign(g)) < 0) absg(g);
        itog(1,scratch3);
        gshiftleft(n,scratch3);
        itog(1,scratch4);
        subg(scratch4,scratch3);
        if(gcompg(g,scratch3) >= 0) subg(scratch3,g);
        if (s < 0) {
            negg(g);
            addg(scratch3,g);
        }
        if (the_sign < 0) {
            negg(g);
            addg(scratch3,g);
        }
    }

    RETURN_GIANT(scratch3);
    RETURN_GIANT(scratch4);
}



int
mersenneinvg(
    int        q,
    giant     x
)
{
    int        k;

    if (!u0)
        u0 = newgiantbits(1);
    if (!u1)
        u1 = newgiantbits(1);
    if (!v1)
        v1 = newgiantbits(1);

    itog(1, u0);
    itog(0, u1);
    itog(1, v1);
    gshiftleft(q, v1);
    subg(u0, v1);
    mersennemod(q, x);
    while (1)
    {
        k = -1;
        if (isZero(x))
        {
            gtog(v1, x);
            return(0);
        }
        while (bitval(x, ++k) == 0)
            ;

        gshiftright(k, x);
        if (k)
        {
            gshiftleft(q-k, u0);
            mersennemod(q, u0);
        }
        if (isone(x))
            break;
        addg(u1, u0);
        mersennemod(q, u0);
        negg(u1);
        addg(u0, u1);
        mersennemod(q, u1);
        if (!gcompg(v1,x))
            return(0);
        addg(v1, x);
        negg(v1);
        addg(x, v1);
        mersennemod(q, v1);
    }
    gtog(u0, x);
    mersennemod(q,x);
    return(1);
}

void
powermod(
    giant        x,
    int         n,
    giant         g
)
/* x becomes x^n (mod g). */
{
    giant scratch2 = BORROW_GIANT();
    gtog(x, scratch2);
    itog(1, x);
    while (n)
    {
        if (n & 1)
        {
            mulg(scratch2, x);
            modg(g, x);
        }
        n >>= 1;
        if (n)
        {
            squareg(scratch2);
            modg(g, scratch2);
        }
    }

    RETURN_GIANT(scratch2);
}


void
powermodg(
    giant        x,
    giant        n,
    giant        g
)
/* x becomes x^n (mod g). */
{
    int         len, pos;
    giant        scratch2 = BORROW_GIANT();

    gtog(x, scratch2);
    itog(1, x);
    len = bitlen(n);
    pos = 0;
    while (1)
    {
        if (bitval(n, pos++))
        {
            mulg(scratch2, x);
            modg(g, x);
        }
        if (pos>=len)
            break;
        squareg(scratch2);
        modg(g, scratch2);
    }
    RETURN_GIANT(scratch2);
}


void
fermatpowermod(
    giant     x,
    int        n,
    int        q
)
/* x becomes x^n (mod 2^q+1). */
{
    giant scratch2 = BORROW_GIANT();

    gtog(x, scratch2);
    itog(1, x);
    while (n)
    {
        if (n & 1)
        {
            mulg(scratch2, x);
            fermatmod(q, x);
        }
        n >>= 1;
        if (n)
        {
            squareg(scratch2);
            fermatmod(q, scratch2);
        }
    }
    RETURN_GIANT(scratch2);
}


void
fermatpowermodg(
    giant     x,
    giant    n,
    int        q
)
/* x becomes x^n (mod 2^q+1). */
{
    int        len, pos;
    giant    scratch2 = BORROW_GIANT();

    gtog(x, scratch2);
    itog(1, x);
    len = bitlen(n);
    pos = 0;
    while (1)
    {
        if (bitval(n, pos++))
        {
            mulg(scratch2, x);
            fermatmod(q, x);
        }
        if (pos>=len)
            break;
        squareg(scratch2);
        fermatmod(q, scratch2);
    }
    RETURN_GIANT(scratch2);
}


void
mersennepowermod(
    giant     x,
    int        n,
    int        q
)
/* x becomes x^n (mod 2^q-1). */
{
    giant scratch2 = BORROW_GIANT();

    gtog(x, scratch2);
    itog(1, x);
    while (n)
    {
        if (n & 1)
        {
            mulg(scratch2, x);
            mersennemod(q, x);
        }
        n >>= 1;
        if (n)
        {
            squareg(scratch2);
            mersennemod(q, scratch2);
        }
    }
    RETURN_GIANT(scratch2);
}


void
mersennepowermodg(
    giant     x,
    giant    n,
    int        q
)
/* x becomes x^n (mod 2^q-1). */
{
    int        len, pos;
    giant    scratch2 = BORROW_GIANT();

    gtog(x, scratch2);
    itog(1, x);
    len = bitlen(n);
    pos = 0;
    while (1)
    {
        if (bitval(n, pos++))
        {
            mulg(scratch2, x);
            mersennemod(q, x);
        }
        if (pos>=len)
            break;
        squareg(scratch2);
        mersennemod(q, scratch2);
    }
    RETURN_GIANT(scratch2);
}



void
fermatmod (
    int             n,
    giant             g
)
/* g := g (mod 2^n + 1). */
{
    int the_sign, s;
    giant scratch3 = BORROW_GIANT();

    if ((the_sign = gsign(g)) < 0) absg(g);
    while (bitlen(g) > n) {
        gtog(g,scratch3);
        gshiftright(n,scratch3);
        subg(scratch3,g);
        gshiftleft(n,scratch3);
        subg(scratch3,g);
    }
    if(!isZero(g)) {
        if ((s = gsign(g)) < 0) absg(g);
        itog(1,scratch3);
        gshiftleft(n,scratch3);
        iaddg(1,scratch3);
        if(gcompg(g,scratch3) >= 0) subg(scratch3,g);
        if (s < 0) {
            negg(g);
            addg(scratch3,g);
        }
        if (the_sign < 0) {
            negg(g);
        }
    }

    RETURN_GIANT(scratch3);
}


#if 1
void
smulg(
    unsigned long    i,
    giant             g
)
{
    if (!i) {
        g->giantDigits = 0;
        return;
    }

    if (!isZero(g)) {
        unsigned long topresult = MultVecsByULong(    g->vectors,
                                                    g->giantDigits,
                                                    i);

        if (topresult) {
            AddMSDLong(g, topresult);
        }
    }

    ASSERT_GVALID(g);
}
#else
void
smulg(
    unsigned short    i,
    giant             g
)
/* g becomes g * i. */
{
    giant    temp = BORROW_GIANT();

    itog(i, temp);

    mulg(temp, g);

    RETURN_GIANT(temp);

}
#endif

long            icompg(unsigned long i, giant g)
{
    long    result;

    giant    temp = BORROW_GIANT();

    itog(i, temp);
    result = gcompg(temp, g);

    RETURN_GIANT(temp);

    return result;
}


void         trunctobitlen(long len, giant g)
{
    unsigned long     neededDigits;
    unsigned long     trimbits;
    unsigned char    *pDigit;
    unsigned char     maskByte;
    long            i;

    if (bitlen(g) < len) return;

    neededDigits = (len + BITS_PER_VECTOR - 1)/BITS_PER_VECTOR;

    g->giantDigits = neededDigits;

    trimbits = len & 0x7f;

    trimbits = BITS_PER_VECTOR - trimbits;

    if (trimbits) {
        pDigit = (unsigned char*)&g->vectors[neededDigits-1];

        for (i=0; i<trimbits >> 3; i++) {
            pDigit[i] = 0;
        }

        trimbits &= 0x07;
        maskByte = 0xff;

        maskByte = maskByte >> trimbits;

        pDigit[i] &= maskByte;

    }
}


#if !INLINE_GIANTS

/* Returns 1, 0, -1 as a>b, a=b, a<b. */
long gcompg(giant a, giant b)
{
    ASSERT_GVALID(a);
    ASSERT_GVALID(b);

    if (a->giantSign > b->giantSign) return 1;
    if (a->giantSign < b->giantSign) return -1;

    return a->giantSign * CompareMagnitude(a, b);
}

signed long        gtoi (giant g)
{
    GAssert(bitlen(g) <= 32);

    if (!g->giantSign) return 0;

    return g->giantSign * ((unsigned long *)g->vectors)[LONGS_PER_VECTOR-1];
}

/* Returns the giantSign of g: -1, 0, 1. */
long             gsign(giant g)
{
    ASSERT_GVALID(g);

    if (!g->giantDigits) return 0;

    return g->giantSign;
}

void         negg(giant g)
{
    g->giantSign = -g->giantSign;
}

void        setsign(giant g, long sign)
{
    if (sign > 0) {
        g->giantSign = 1;
    } else {
        g->giantSign = -1;
    }
}

void        absg(giant g)
{
    if (g->giantSign < 0) g->giantSign = 1;
}

long isZero(giant g)
{
    ASSERT_GVALID(g);

    return (g->giantDigits == 0);
}

Boolean    isone(giant g)
{
    Boolean                    result = false;
    unsigned long             *pLongPtr;
    unsigned long            accumulator;
    long                    counter;

    if (g->giantSign == 1 && g->giantDigits == 1) {
        pLongPtr = (unsigned long*)g->vectors;
        accumulator = 0;
        for (counter=0; counter < LONGS_PER_VECTOR -1 ; counter++) {
            accumulator |= pLongPtr[counter];
        }

        // if first three longs are zero
        if (!accumulator) {
            // and last one is one, then it's one
            result = (pLongPtr[counter] == 1);
        }
    }

    return result;

}

Boolean            isOdd(giant g)
{
    return (GetNthByte(g, 0) & 0x01);
}


#endif //#if INLINE_GIANTS